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1 . 


Introduction 


One of the principal objectives in the future utilization of 
orbital laboratories and stations is the production of materials 
in space. This type of material production is referred to as space 
processing. One of the primary constituents of most production processes 
involves the solidification of materials. The broad objective of the 
present report is to investigate the solidification process in a low 
gravity environment, while the specific task is to look at the stability 
of the solid-liquid interface under the same conditions. The meaning 
of the term stability as used here will be discussed more fully later 
on in the report. 

The problem of material solidification in general has attracted 
researchers for many decades due to both its importance to industry 
from the practical point of view and its mathematical difficulty from 
the theoretical point of view. The analytical solidification problem 
is commonly known as the "Stefan Problem" in honor of Jakob Stefan 
(1889), who first formulated it. At the present time, solutions to 
various aspects of the problem exist under different boundary condi- 
tions, such as to make it of practical use for industrial applications. 
For recent reviews of such solutions the reader is referred to the 
books by Rubinstein (1971) and Ockendon and Hodgkins (1975). Also, 
the not-so-recent review articles by Boley (1969) and Muehlbauer and 
Sunderland (1965) are informative. However, the problem in its entirety 
is far from solved, as is revealed from a scan of the recent literature 
which, incidentally, also shows the continued interest in the problem. 


The mathematical difficulty of the problem is manifested in the 
fact that although the governing equations are linear, the boundary 
conditions are strongly non-linear, with the added difficulty that 
the boundaries are not stationary. Due to this second complication, 
almost all of the solutions to the problem are time dependent with 
only very few exceptions where it is steady. 

Careful inspection of the problem reveals that the formulation 
of the problem is not at all altered when the surrounding environment 
is at reduced gravity. Hence, all of the mathematical difficulties 
alluded to earlier are present under the new conditions. However, 
instead of tackling the full-blown problem for the present report 
and thus adding relatively little new information, it was decided 
to look at one aspect of the problem which might be of relevance to 
space processing. This aspect involves the deformation of the solid- 
liquid interface from its original form during a solidification process 
in a reduced gravity environment. 

The problem just defined above is not unique to space process- 
ing, and there exist quite a few applications for ground base work. 

If the solidifying material is composed of a binary alloy, for example, 
it has been shown first by Mullins and Sekerka (1964) and later by 
Wolkind and Segel (1970) that the solid-liquid interface could deviate 
from its original planar configuration for both finite and infinitesimal 
disturbances. On the other hand, if the solidification process is 
allowed to proceed from the upper boundary of the liquid, instabilities 
of the Benard type might contribute to the deformation of the 
phase boundary. Such a problem has been investigated by many workers, 
some of whom are Turcotte (1974), Schubert et al . (1975), and Busse 
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and Schubert (1971), where the application was to mantle convection. 
Also, the same type of instabilities might set in during a melting 
process if the melting is from below, as is shown analytically by 
Sparrow et al . (1977) and later verified experimentally also by Sparrow 
et al . (1978). 

In the present report the stability of the solid-liquid interface 
of a pure metal will be studied during a solidification process in 
which the liquid surface is free. Two such problems will be considered, 
one in which the phase boundary is stationary and another in which it 
is propagating at a constant speed. Note that the added complication 
due to the solidification of a mixture or a binary alloy is avoided 
here since it will not greatly enhance the basic understanding of the 
problem. Since the problem will be analyzed in a zero gravity environ- 
ment, buoyancy effects do not contribute to the problem, thus leaving 
only the influence of a variable surface tension to be investigated. 

Seki et al . (1977) have studied the influence of surface tension 
coupled with buoyancy effects on the stability of the phase boundary. 
Their work, however, relates to a stationary phase boundary. Chang 
and Wilcox (1976), on the other hand, considered the effects of surface 
tension in a floating zone melting process. However, their analysis 
was not of the stability type but rather numerical in which the 
solidification process does not enter. 

In the next few sections of this report the problem to be studied 
will be defined and the equations will be set up with the appropriate 
boundary conditions. Then the problem will be solved and the results 
tabulated and discussed. In the final section these results will be 
discussed as they apply to real working materials such as some pure 
metals whose use is anticipated in space processing applications. 
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2 . 


Mathemati cal Formul ati on 


The complete analysis of the shape of the solid-liquid interface 
in a general solidification process is extremely complicated, which 
ultimately might require a complicated numerical solution and thus 
is beyond the scope of the present report. Alternatively, it is 
possible to study the stability of the interface for a simple geo- 
metrical configuration. By stability of the interface the question 
is asked: given an original shape of the interface, say, planar, or 
any other configuration, under what conditions will the shape remain 
unaltered during the various stages of the solidification process? 

This is a typical example of a classical hydrodynamic stability problem 
in which it is determined whether the interface geometry will change 
if a slight disturbance is imposed on the original configuration. 

The evolution of this disturbance with time will then determine the 
fate of the original configuration. If the disturbance decays with 
time, then the original confi guration is said to be stable with respect 
to the disturbance, and it is unstable otherwise. If the disturbance 
is infinitesimal, then the stability analysis is called linear, and if 
it is finite then the analysis is nonlinear. In the present report, 
only linear stability analysis will be performed. 

Since most stability analyses are performed with respect to time, 
it is imperative that an original state is found which is either 
independent of, or at most periodic with, time. As discussed in the 
introduction, in a solidification process the solution is almost 
always a nonperiodic function of time. For the present study two 
solutions for the mean state were found which are independent of time. 
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both of which are analyzed here. The first configuration which is shown 
in Figure 1 constitutes a fluid and a solid layer which are in equilib- 
rium and hence the liquid-solid interface is planar and stationary. 
Although such a configuration is not strictly a solidification process 
it is nevertheless used for illustrative purposes and further to 
unravel the main difficulties associated with the problem. The second 
configuration constitutes a continuous solidification process in which 
the interface boundary is moving at a constant rate. With a simple 
transformation it will be shown that a steady state solution exists 
for this process, thus facilitating a managable linear stability 
analysis. This second example is closer to a real solidification 
process and thus is a more realistic example than the first. 

Since the present stability analyses are linear, a direct answer 
to the stability of the interface cannot be obtained. The stability 
or the instability of the interface will be considered as a conse- 
quence of the stability or instability of the fluid. Thus, if the 
analysis reveals the liquid to be unstable to infinitesimal distur- 
bances, it is then assumed that convection of some sort will set in 
the liquid in such a way as to disturb the original shape of the 
interface and probably deform it, as shown schematically in Figure 2. 

If the fluid is found to be stable, on the other hand, it will be 
assumed that any initial deformation of the interface will be evened 
out and the interface will retain its original configuration. 

2.1 Steady State Solidification Process 

A. The Mean State Solution 

It can be formally proven, and will be shown below, that if the 
upper surface of a melt layer which is on top of a solid layer whose 
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surfaces are kept at a constant but different temperature each, then 
part of the melt will solidify until the solid-liquid interface reaches 
an equilibrium position. This equilibrium position of the solid- 
liquid interface is a function of both temperatures and the solidifica- 
tion temperature only. This is one of the basic and simple steady 
state solutions to the Stefan problem [Rubinstein (1971)]. As a first 
example, the stability of this solution will be studied where a sketch 
of the problem is provided in Figure 1. 

Consider two finite depth layers, one of melt and another of 
solid, which are infinite in the horizontal direction, where the 
temperature distribution in both the liquid and solid is denoted by 
T^ and T^, respectively. Since there is no convection in the mean 
state, the energy equations for a steady state are given by 

d ^ os ^ y d “Kt ^ O J 

dV'' ( 1 ) 

with the temperature at the top of the liquid layer given by 
To£(d)=T 2 and the temperature at the bottom of the solid is given by 
Tog(0)=T^. In here the asterisk denotes dimensional quantities. 

However, since we are considering a solidification process, the energy 
balance and conservation of mass considerations at the solid-liquid 
interface imposes the following three additional conditions: 

> r 

4s * 

/: d 7-s r k d'TJi 
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where k and k. are the thermal condictivi ties of the solid and liquid, 
s ^ 

★ 

respectively. is the solidification temperature and S is the 
equilibrium position of the interface. 

A solution to equation (1) with all the boundary conditions is 
given by 

(3) 

i 


The equilibrium position of the interface is uniquely determined by 
this solution to be 

■T; - KTi ^ 

where (4) 

/s = ^ h 

The solution given by (3) is identical to that obtained by Rubinstein 
(1971) for the time dependent problem in the limit as time->-“ . Further- 
more, the above solution is no more than the steady conduction solution 
for the fluid and solid layers with a linear temperature distribution 
in each. 
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B . Linear Stability Analysis 

In order to investigate whether the interface boundary can withstand 
a slight perturbation, a linear stability analysis was performed on 
the mean state solutions (3) and the boundary conditions. The analysis 
was carried out along lines which are similar to what is commonly done 
in hydrodynamic stability analyses. The field variables are first 
decomposed into a mean and a perturbation in the form 

7 ^ - ■«' / > 

ZT ^ ^ u' ^ (5) 

where T is the temperature field and u is the velocity vector, in 
the present case is zero since it is assumed that there is no convective 
motion in the mean state. These functions are then substituted into 
the momentum and energy equations and the mean state is subtracted 
out. Since we are considering only a linear stability analysis, the 
remaining equations for the perturbation functions are then linearized. 
Furthermore, these equations are then manipulated to eliminate the 
pressure term resulting in the following equations; 
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Use has been made in the above equations of the mass conservation 
equation. In these equations v is the kinematic viscosity of the 
liquid, K and k are the thermal di ffusi vities of the liquid and 

X/ w 

solid, respectively, and w' is the z component of the velocity vector, 
d T 

and o£ is the temperature gradient of the mean state in the liquid 
dz 

which is a constant in this case. The independent variables in equa- 
tions (6)-(8) have been made dimensionless in the following way: 

( ^ /d . 1 . 


and 



' 'd /^i 


V 







( 9 ) 


where again the asterisk indicates dimensional quantities. 

The set (6)-(8) constitutes a linear partial differential system 
for the perturbation variables w’ , T^ and T^. A solution to this set 
is sought for appropriate boundary conditions which are obtained in the 
following way. At the lower surface of the solid it is assumed that 
the temperature perturbation is zero, i.e.. 


r^"'(o) = o 

At the top of the surface of the liquid we require the velocity 
perturbation to be also zero, i.e.. 


(10) 


ur * ' { j + ^ ) s o 


(11) 


9 


Also, a radiation condition is imposed on the upper surface of the 
liquid in the form 


where q represents the rate of change with temperature of the rate of 
loss of heat per unit area from the upper surface to its upper environ- 
ment. Note if q^co we approach a conducting boundary condition, while 
if q->o we approach an insulating boundary condition. A comprehensive 
discussion of this boundary condition is provided by Pearson (1958). 
Another condition to be imposed on the upper surface of the liquid 
layer is that of the balance of the tangential forces on that surface, 
which is given by Levich and Krylov (1969): 


(12) 



where a and p are the surface tension and dynamic viscosity of the fluid. 
Now, since the surface tension is known to be a strong function 
of temperature, the above condition will be valid only if the tempera- 
ture of the upper surface is not constant. If instabilities will set 
in the fluid, the temperature distribution will not be constant on the 
surface and the surface tension force will play a major role in the 
stability analysis. For the purpose of the present report it will be 
assumed that the surface tension is a linear function of temperature 
in the form 
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where b = and a is a constant. With this approximation the 

condition foi^' the balance of the tangential forces may be rewritten as: 



There remains the boundary conditions at the interface to be 
satisfied. The first is the requirement that the temperature at the 
interface regardless of the position of the interface remains the 
solidification temperature leading to 


Here it is assumed that the perturbed position of the interface to a 
first order of approximation is the same as its unperturbed position. 
Also, at this position, where the solid boundary acts as a wall, the 
nO“Sl i p conditions are imposed and take the form 

- o ^ ^ ^ , 

35 * 

while the linearized form of the energy balance at the interface leads 
to 


(14) 


(15) 


' _ fe. 

where p is the density of the solid and L is the latent heat of fusion, 
s 

Condition (16) is not intuitively obvious and is rigorously derived 
in Appendix 1 for the present case and for the subsequent case also. 

The linear system of equations (6)-(7) together with the separable 
linear boundary conditions (10)-(16) admit solutions through the 



= - ^ A 


(3 t S 


11 



separation of variables method. To obtain such a solution the perturba- 
tion variables must take the following functional form: 

[ W" ' Tj"' . 7^ ] = 9^^i) X . 

(17) 

Note that only a two-dimensional form of the perturbation function is 
considered, where it is hoped that the analysis is not affected in a 
major way. 

Substituting (17) into (6)-(7) leads to the following system of 
ordinary differential equations for the perturbation functions 


[q, _ (18) 


[ a. - A = 


(19) 


where 



• ^or + ’ i 

The question of stability or instability can be settled by deter- 
mining the value of a^. Specifically, only the sign of ReCa^) is needed 
since if Re(a^) is negative, zero or positive then the system is stable, 
neutral or unstable, respecti velly. To even determine the sign of Re(a^) 
is quite tedious and instead the principle of exchange of stability is 
invoked here. This principle implies that if Re(a^)==0 then lm(a^)=0 
and if this is true, then it is possible to set a^=0 in order to deter- 
mine the neutral stability criteria. Adopting this principle reduces 
the system of equations to the following: 
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( 20 ) 


{ b* - 6:4 ^ O 

while the boundary conditions will take the following form: 


(') ♦ 6 (9| (1) 5 o 


c? 

(o) s <0 
(O) ^ O 


- ^* ) iO) c. O 

')> ^,(0) s K (o) - ( c>’_ )<^ (°) 


(9s (~/i ) ^ o 




where 


( 21 ) 


( 22 ) 


i?> -- f J/t, 

Z'/ — ( / 3T ){ '^lot d 


<^ps (Ti-t:) 

In the above boundary conditions the origin has been shifted to the 
interface boundary. 

The solutions to equations (.20) and (21) may be written immedi- 
ately in the form 


^ a) f /?, + -t ) Co<Jk C» i 

+ ( 4 I A ij 4 -l‘ ^ s.^J^JL, oc I 


(.23) 
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( 24 ) 




r 4 cx^ H C< ^ 


Satisfying the boundary conditions for the perturbation temperature 
in the solid leads to 

and thus 

<pj = o . 

However, A-| through Ag may be determined from the conditions on the 
perturbation temperature in the liquid. Since we are looking for 
neutral stability curves the exact values of A^ need not be determined 
but if all the conditions are to be satisfied and since the conditions 
are homogenous, the requirement that the coefficients be unique will 
lead to the following parameter relation: 


(25) 




X ( SjitJL < ) oc 


(26) 


where , D2 and are the following determinants: 




.SlcU Ok 

' J 




zsa:JL k 
- £>< 






ZpK tX^fCpi 


5 \<i< Ca^K 


( 27 ) 
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2. o< . K 

II 

o 




^ Lf 


K OC 

-z£^ -Km. 




rr 

A Y 


(28) 






o 


p< SinJi^O{ ) +1 K 






1^ K tx k 

This equation gives the requiretd relationship between the Marangoni 
number and the wave number of the neutrally stable disturbances. 

2.2 Constant Rate Solidification Process 
A. The Mean State Solutions 
The analysis performed in the previous section was for a stationary 
interface in which the solidification process entered only in the per- 
turbation equations and not in the mean solution. Thus, the example 
does not realistically model a solidification process and was used 
only for the simplicity of the mathematics involved and for illustra- 
tive purposes. In this section the analysis of a more realistic 


AY J 

(29) 
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problem is presented. The problem will involve a continuous solidifica- 
tion process in which the melt is being fed continuously and at a 
constant rate while the solidification rate is proceeding also at a 
constant rate. A sketch of the problem is shown in Figure 3. However, 
in order to keep the solidification rate constant for a continuous 
process, it is required that the heat removal rate at the bottom of 
the solid be also constant. The physical analog of this model could 
be a continuous solidification process in which the melt is being 
replenished all the time in an exact amount to keep its mass constant 
for a very slow solidification process. 

Thus, consider a finite depth melt which is infinite in the hori- 
zontal extent on top of a solid layer. The temperature distributions 
in both the solid and the liquid are given by the one-dimensional energy 
equations [Carslaw and Jaeger (1959)] in the form 



where p and p are the densities of the solid and liquid, respectively, 
and the rest of the notation is the same as that used in the previous 
section. Note that here the interface position $(t) is a function 
of time and not a constant as before. However, it will be assumed from 
the start that the interface moves at a constant speed v^, where 

dS - cr 
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At this point, as will prove beneficial later, the origin of the 
coordinates is moved to the solid-liquid interface via the following 
Gallelain transformation: 


t 



i 


K 







This transformation will render equations (30) and (31) to the 
following form: 


( 33 ) 



dl” 









(34) 


(35) 


Also in these equations the primes have been dropped for ease of notation. 
The boundary conditions appropriate to this problem are the following: 


'os “ 

d7^5 - It- J - 





(36) 


■^1 





t = d 


(37) 




( 38 ) 
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Note that for the solid boundary condition a heat transfer rate condition 
is used rather than a fixed temperature condition in order to be 
consistent with the assumption of a constant solidification rate. 

The temperature distribution in the liquid phase may be obtained 
from the solution of the above system to be: 


Zf = r * - )][ J - ei*p 


_ / 


where and z are dimensionless quantities having the following form: 


( 39 ) 




— 9 - f 

z u, 




A plot of this distribution of temperature is given in Figure 4. The 
system of equations and the boundary conditions also specify uniquely 
the temperature distribution in the solid and the propagation speed 
of the interface. However, use will be made of the latter condition 
which is also a matching condition and will be written out when the 
need arises. 

B. Linear Stability Analysis 

As indicated earlier, the stability of the solid-liquid interface 
is inferred from the stability of the liquid phase. The stability 
analysis here is performed in a very similar manner as was done in 
the previous section. Again, the temperature and velocity fields are 
separated into a mean and a perturbation and substituted into the 
momentum and energy equations resulting in a linearized form for the 
equations for the perturbation functions. In this case, however. 
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a Gallehian transformation is again adopted as was done for the mean 
temperature distribution where 



Also, the resulting independent variables in the new frame of reference 
are non-dimensional i zed in the following way: 

t -- tV,./j 

L 3 = L x/j ^ 

u ‘ u “/ir 

' r 

where d is the depth of the liquid layer. 

The resulting linear non-dimensional governing equations for the 
perturbation functions are 



Most boundary conditions under which the above equations are to be 
solved in the moving frame of reference are identical to conditions (10)- 
(16) except for condition (16) which now is given by 
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I I Hill 


IIIIHI II 


■ Mill 


III! 


k 


- k, ?Je 

' Pi 


-^^Lur 


Again, due to the form .of the boundary conditions , equations (40)- 
(42) may be solved through the normal mode method. It will be assumed 
for simplicity that the perturbation functions are two-dimensional. 
Under these circumstances the perturbation functions will take the 
same form as (17). Substituting this functional form of the solution 
into equations (40)-(42) we get after setting = o: 



J. 

/? 






o 


(43) 


.r («1 


J 

where 

(? = 

In equations (43)-(45) we have again invoked the principle of exchange 
of stability. 

The solution method for the system of equations (43)-(45) is 
slightly different here from that of the previous case. First, 
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equation (43) is solved alone for the velocity perturbation v(z). 
Then equation (44) is solved for the temperature perturbation 9j^:(z) 
as an inhomogenous equation. Equation (45), however, is uncoupled 
and may be solved separately. The solution for equation (43) may be 


easily written in the form 


/I; 




where are given by 


(46) 


'X.,^ = t 


X 


while the solution to equation (44) is given by 

uc 

where through and X-j through are the same constants as those 
in the velocity distribution, while Ag and Xg are given by 

The Cj appearing in (48) are given by 


(47) 


(48) 


(49) 




i Si j 4 


( 50 ) 
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where 


^ { / - e^pi- 

>>«, 

Since the solution for the perturbation velocity and temperature are 
obtained separately it is instructive to give the boundary conditions 
in terms of both variables. These are 
^(O) iC? 

4 fo) = O 
fo) r: O 

^ ^ (f (o) ^ o 

K 

(rd) ^ o i ^ 

r> Of ti) ^ - o 

)(T it] (i) ^ a 

where 


^ = /\pe.CK 

^ z ( ) 

C 

Thus substituting the solutions given by (46) and (48) into the boundary 
conditions (52) yields a system of equations for A. If the value of 
the coefficients is to be unique then the determinant of the coeffi- 
cient matrix for A^. must be zero, i.e.. 


L£ ] = O 


where the matrix E is given by 



C. 


C c 

i / 


■^t 


0 0 


e/jrp( 

€>ypf ^3) €/yp{^ii^ 

0 0 

i 



^S- ^C 

€./jip 

6-^p ^'Xk»+^i1 

1 '^Kl'+ ^ i. *^w>'f '^\|H 

■^f t^Hl‘*'^\) ^P l.^Vk+'^«l^ 


Z« e^pt ^, ) 

C^p*X^^y^X (^j C ^l4 ^ ( ^s\ 

^p\ 

— 



J 


( 55 ) 



To obtain the neutral stability curves equation (54) is then solved for 
one of the parameters in terms of the rest, i.e.. 




( 56 ) 
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3 • Non-Dimen sional Parameters 

The goal of the present work is to calculate the neutral 
stability curves and the critical Marangoni number, M, as a 
function of the non-dimensional wave number of the disturbance. 

It was necessary to non-dimensionalize the problems to be able 
to make these calculations. There are always many ways to 
non-dimensionalize a problem, but the method used recover 
familiar dimensionless parameters plus additional parameters. 

3.1. Steady State Process 

The dimensionless independent variables were defined along 
with the solution to the problem given earlier. The dependent 
variables were non- dimens ionalized to give the following relation- 
ship for the neutral stability curve: 


M = M(a,B,A,Y,A, K) 


where 


* 

T^>r _ da *^o£ j2 
^ - dT* di^ ^ 




2-Kd 


a = dk = — ^ = dimensionless V7ave n\amber (X = disturbance 
wave length, k = disturbance wave number) 


B = 3^ = Biot modulus 


A = 




Cp^(VTi) 


dT 


Y = 


oil 


dz 
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It was desirable to calculate the neutral stability curves 
using values of these parameters which characterize material 
systems of interest for space processing. Table I lists property 
values for some metallic materials that would be of interest. 

The parameters A and K depend only upon the material properties 
and their values are tabulated in Table II. A varies from 0.7 
to 3.0 for the metals considered, while K varies from 0.44 
to 1.74 but. is generally around 0.5. 

The remaining dimensionless parameters depend upon the 
boundary conditions of the problem as well as the material 
properties. First consider the mean temperature field. The 
analysis assumed that the material properties were constant 
but in fact they are strong functions of temperature. There- 
fore, for this reason and because solidification usually 
proceeds slowly, the temperature difference across the liquid- 
solid system can be assumed to be small, say, 2®C. Then A 
varies from 35 to 200 and y from 1.2 to 1.5. The Biot modules, 
B, on the other hand, which characterizes the heat transfer 
through the upper liquid surface, can be left as a parameter 
in the problem. 


^Smaller values of A correspond to large values of (T*^-t|) and 
cannot be expected to yeild satisfactory results using the 
present analysis. 


3.2. Con stant Rate Process 

For the continuous rate process the neutral stability 
curve is given by the following relation: 

M = ii^(a,R,Pr,Q,B,A,K,P,C) 


where 


R = 




= Reynolds number 


V 


Pr = 


- Prandtl number 


I” = 


C = 


Q = 


^qsd 

^ ^ -k 
(T«-T ) 
^ 2 m 


A = 


^ 7C /V . 

Cp (T2-V 

s 


and K and B are defined as before. Notice that A is now 

defined in terms of T^-T* and is thus different from the parameter 

2 m 

defined for the steady state problem. M is defined as 


* da 


M = 


(To-T^) 


2 dT’ 


If a mean temperature gradient in the liquid defined by 




dz'" , d 
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is used in the definition of the Marangoni number, then M 
is related to it as follows; 


M = R Pr M 


However, M is the more natural parameter for the analysis and 
the neutral stability curves will be given in terms of it. 
Since Pr is small M is considerably smaller than M, for 
Reynolds numbers of the order of one. 

The continuous rate problem was non-dimensionalized using 
the* speed of the phase plane, Vp , as one of the parameters in 
(M and R) . The phase plane speed is not one of the boundary 
conditions of the problem but rather is part of the solution. 
It is determined from the non-linear matching condition that 
results from matching the solutions at the liquid-solid 
interface for the mean motion. That matching condition is the 
following : 

(tZ-TJ) 








= PsLUp. 


This relation can be non-dimensionalized using the parameters 
listed above and they were selected, in part, to accomplish that 
non-dimensionalization . In non-dimensional form the equation 
becomes 


QPC 

ARP*- 


ex p 


(-¥^) 




A 1^1 - p f-R P Pr)] 


-2 


From Table II it can be seen that Pr = p(10 ), A=0(1) , and 

P = 0(1). Therefore, the exponentials are small if R = 0(1). 

-2 

R = v„d/v„. From Table I, v. = 0(10 ) and for a typical 
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_ 2 

system d = 0(1 cm) and Vp will be of the order of 10 cm/sec if 
R is of the order of one. 

Therefore, for R = 0(1) the above expression can be 
linearized, yielding the following expression for the Reynolds 
number in terms of material properties and the boundary condi- 
tions of the problem. 

R = PC I Q+1 V 
^ Pr V QK+A ) 

Values of the parameters for some metallic materials are given 
in Table II. The above expression is used to calculate R in 
terms of the other given parameters. Pr, P and C depend only 
upon the properties of the materials. Q, which is the non- 
dimensional heating rate, fixes the rate of solidification and 
can be independently varied. B, which determines the perturbation 
heat transfer rate, can also be kept as a parameter in the 
solution. 
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4. Discussion 


It is a common practice in stability analyses to present 
the results in the form of neutral stability curves such as 
shown in Figure 5. Such curves delineate in the relevant 
parameter space regions of stability and instability. It 
should be emphasized that the information obtained from such 
curves is the most that one can hope to get from linear stability 
analyses. However, it should be remembered that in obtaining 
such information the principle of exchange of stability has 
been used without any proof on its validity for the present 
problem, although it is felt that it may hold here. For more 
information on the ensuing process, whether in the form of a 
stationary time periodic bifurcated solution or a totally 
unstable solution, resort must be made to nonlinear analysis. 
Thus, the information obtained here is just qualitative in 
nature without any detail on the ensuing flow, if it occurs. 

4. 1 Steady State Process 

The neutral stability curves for the steady state solidi- 
fication process given by equation (26) are shown in Figure 5. 

All of these curves show a critical Marangoni number indicating 
a range of Marangoni numbers, M, for which the perturbations 
are stable regardless of the wave number. The various curves 
are of different values of the Biot modulus, B, which is a non- 
dimensional form of the rate of heat transfer through the upper 
liquid surface. The similarity between this figure and Figure 1 
of Pearson (1958) is striking. However, this similarity is 


30 


expected since the problem solved here is almost identical 
with that of Pearson. The only difference between the two is 
in the form of the energy balance condition (16) , whereby 
Pearson's condition is recovered when the latent heat parameter 

A . In fact, this is exactly what happens, for a's A becomes 

large, Pearson's stability curves for the ideally conducting 
case are recovered. However, this result is achieved quite 
fast and there exists differences only for unrealistically 
small values of A, i.e. , A values of order 1. For any realistic 
solidification process in which the temperature gradient in the 
liquid phase is not too high, the value of A is of order 10 to 
100. The conclusion to be drawn here is that the system becomes 
more unstable (shifting of the neutral stability curves to the 
left) with a decrease in the latent heat values. 

4 . 2 Constant Rate Process 

The neutral stability curve is given as the function 

M = M(a) 

with Pr, Q, B, A, K, P, and C as parameters. Pr, K, P, and C 
depend only upon the properties of the materials and neutral 
stability curves were computed for typical values of these 
quantities (see Table II) . The dependence of the neutral sta- 
bility curve on each of these parameters was examined and will 
be described below. Although the problem examined is greatly 
simplified compared to a realistic time- dependent solidification 
process, it is hoped that the dependence of the stability upon 
these parameters will be correctly predicted. All results are 
tabulated in Appendix III. 
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Increasing Q by increasing qos » which increases the rate 

of solidification, leads to lower critical values of M. Typical 

neutral stability curves are shown in Figure 7. The values of Q 

used correspond to very small values of ^qs * order of 0.1 

2 

to 0.5 joules /sec cm . It is of interest to note that the wave 
number corresponding to ^^^rit approximately constant (except 
as noted later) and equal to 2.0 to 2.2. This corresponds to 
a most critical wave length of about ird. This can be expected 
to be the dimension of the cellular pattern that would occur 
when ^crit exceeded (Pearson (1958)). 

Increasing the rate of heat transfer through the upper 
liquid surface, i.e., increasing B, stabilizes the flow (Figure 
8) . This conclusion is identical to that found in the steady 
state problem and by Pearson (1958) for a non-solidifying liquid 
Furthermore, this latter result has been borne out by experiment 
Notice now, however, that the wave number for the most critical 
disturbance increases with heating rate (Figure 9) and the con- 
vection cell size can be expected to decrease accordingly. B 
can be varied by changing either the depth of the liquid or the 
heat transfer coefficient, q, and its variation has no effect 
on the solidification rate R. 

A can be changed in two ways and they have different in- 
fluences upon the staiblity of the system. Different materials 
will have different values of L/Cp^ (Table II) and this will re- 
sult in different values of A, for a given (TJ-T*) . Keeping all 
other boundary conditions constant still results in a change in 
R. The effect of the change in L/Cp upon is shown in 
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Tabel IV and Figure 10. Materials with higher values of L/C^^ 

^s 

are more stable, even though the rate of solidification is also 
in general faster. 

A can also be varied by changing (T^-T*) but that also 
changes Q. The effect of changing (TJ-T*) in A and Q simul- 
taneously, while keeping the other boundary conditions constant, 
is shown in Table IV and Figure 11. Again, increasing A, through 
decreasing (TJ-T*) , stabilizes the flow. 

As noted above, the parameters Pr, K, P, and C are con- 
stant for a given material. However, the sensitivity of the 
stability upon these parameters was examined so that comparison 
could be made of the relative stability of different materials. 

The Prandtl number can vary over many orders of magnitude, 
depending upon the material (Table II). Because the Reynolds 
number is directly proportional to Pr, calculations were per- 
formed at constant R = 1.0 by varying A simultaneously with Pr. 
The variation of the critical value of M with Pr for this case 
is shown in Figures 12 and 13. Because the Maragoni number is 
related to M by 

M = R Pr 

the variation in critical Marangoni number is very small. The 
results are summarized in Table V. These results are not 
accurate at the larger Pr because very large (T*-T*) would be 
required to make A small enough to keep R = 1.0. 

The effect of density and specific heat ratio, P and C, 
upon the critical value of M is shown in Figures 14 and 15. The 
greater the increase in the density of the material upon 


33 



solidification, i.e., the lower the value of P, the more stable 
the flow. Materials such as bismuth and water with values of P 
greater than one are somewhat less stable than materials with 
P less than one. Note that R increases somewhat as P is in- 
creased verifying the previous conclusion. The same conclusion 
can be drawn for increasing specific heat ratio, C (Figure 15). 

C also influences R but its effect is greater than the effect of 
changing R alone by changing Q. 

K, the ratio of the coefficients of thermal conductivity, 
can be shown to have no effect upon the stability of the system. 
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5.0 Conclusions 


The following conclusions can be drawn from the solution 
of these problems: 

1. Increasing the heat transfer rate from the upper sur- 
face of the liquid by increasing the Biot modulus , B, stabilizes 
the flow. The wave number corresponding to the most tinstable 
disturbance increases with B causing a corresponding decrease 

in convection cell size. 

2. Increasing the rate of solidification by increasing 
the heat transfer rate through the solid (i.e., increasing Q 
by increasing q^g) destabilizes the flow. As Q increases the 
Reynolds number, R, is increased. 

3. The wave number of the most unstable disturbance is 
about 2,0 to 2.2 for B = 1.0 and’ varies appreciably only with 
changes in B. 

4. For a given temperature difference across the liquid, 
(T*-T*) , materials with larger values of L/Cp^ are more stable 
than those with smaller values. 

5. Decreasing (T*-T*) while keeping other boundary con- 
ditions constant stabilizes the flow. 

6. The Prandtl number, Pr, has little effect upon the 
critical Marangoni number for constant R. 

7. Materials with larger values of density and specific 
heat ratios, P and C, are less stable than those with lower 
ratios . 

8. The thermal conductivity ratio, K, has no effect on 
the stability of the flow. 
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Figure 4. The Non-Dimensional Temperature Profile 
in the Liquid for the Constant Rate 
Solidification Process 
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Figure 7, Effect of increasing the rate of solidification on the stability 
of the constant rate process. = 0.02, K = 0.5, P - 0.965, 

C = 1.04, A = 350, B = 1.0. 



W pH 


Figure 8. Effect of increasing the upper heat transfer rate on the 
tability of the constant rate process. = 0.01, K - 
= 0.975, C = 1.01, A = 175, Q = 1.0- 
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Figure 9. Variation of most unstable wave number 
with upper surface heat transfer rate 
for the constant rate process. Pr = 0.01, 
A = 175, Q = 1.0, K = 0.7, P = 1.0, 

C - 1.0, R = 1.14. 





Figure 10. Effect of change in A (with constant (T 2 * 
Tjjj*) upon Afcrit for the constant rate 
process. Pr = 0.01, Q = 2.0, K = 0.7, 

P = 1.0, C = 1.0, B ^ 1.0. 




Figure 11. Effect of change in (T 2 * - (with 

qos constant) upon ^crit the 

constant rate problem. = 0.01, 

K = 0.7, P = 1.0, C = 1.0, B = 1.0. 









p^ 


Figure 13. Variation of A^crit with for R 
for the constant rate process. 

Q = 1.0, K = 0.7, P = 1.0, C = 1 
B-1.0. M = RP^A/. 




Figure 15. Effect of specific heat ratio on Merit 
the constant rate process. 

Pr = 0.02, A = 175, Q = 1.0, K = 0.5, 

P = 1.0, B = 1.0. 



TABLE I. PROPERTIES OF MATERIALS 
a) Properties of Liquid State 


Materi al 


Specific 

heat 

f joul es^ 

( g°c ; 

Dynamic 
viscosity 
(g/cm sec®C) 

i Kinematic 
viscosity 
(cm^/sec) 

Thermal 
conductivity 
/ joules \ 
^sec cm°cj 

Thermal 
di ffusivity 
(cm2 /sec) 

Temperature 
at which 

properties given 
(°C) 

Melting 
1 Temperature 
■ (°C) 

Ref. 

Sodium 

0.928 

1 .378 

6.80x10“^ 

7.34x10“^ 

0.869 

0.680 

97.82 

97.82 

a 

Potassium 

0.828 

0.814 

5.67x10”^ 

6.85x10'^ 

0.541 

0.802 

63.2 

63.2 

a 

Bismuth 

10.09 

0.150 

1 .66x10“^ 

1.64x10'^ 

1 

0.146 

0.096 

304 

271 

b 

Lead 

10.51 

0.147 

2.12x10'^ 

2.01x10'^ 

0.151 

0.098 

400 

327.4 

b 

Tin 

6.98 

0.255 

1.91x10"^ 

2.74x10'^ 

0.304 

0.171 

240 

231.9 

b 

Zinc 

6.92 

0.501 

3.17x10"^ 

4.58x10“^ 

0.577 

0.166 

500 

419.5 

b 

Water 

0.9998 

4.218 

1.79x10"^ ' 

1 

1.79x10"^ 

5.66x10’^ 

1 .32x10“^ 

0.0 

0.0 

c 

























































b) Properties of Solid State 


Material 

Hi 

HwnBSsilH 

Specific 

heat 

1 /joules\ 
: 1 g°c 1 

Thermal 
conductivity 
f joules \ 

1 Isec cm°C| 

i^HI 

Temperature 
at which 

properties given 
{“0 

Latent heat 
(joules/qm) 

Reference 

ilHBIIIIIi 

0.951 

1.361 

1 .191 

0.914 

97.82 

113.2 

a 

Potassium 

0.849 

0.806 

1 .057 

2.45 

63.2 

1 

59.3 

i 

a 

Bismuth 

9.76 

0.125 

0.084 

0.069 

25 

50.2 

c 

Lead 

11.34 

0.142 

0.346 

0.215 

25 

1 

24.6 

c 

Tin 

7.16 

0.226 

0.640 

O'. 396 i 

25 

21.1 

c 

Zinc 

7.40 

0.389 

1 .150 

0.400 

25 

102.1 

c 

Water 

0.917 

2.100 

2.21x10'^ 

-2 

1 .24x10 

0.0 

333.5 

c 


ui 

rsj 

References: a) Foust, 0. J. , editor, Sodium-NaK Engineering Handbook, Volume 1, Sodium Chemistry 

and Physical Properties , Gordon and Breach ^1972). 

b) Ukanawa, Anthony 0., "Diffusion in Liquid Metal Systems," NASA Contractor Report, 
Contract NASS- 30252, June 1975. 

c) Handbook of Chemistry and Physics , 54th edition, CRC Press (1973). 



















































Table II 


VALUES OF NON-DIMENSIONAL PARAMETERS 


Material 

Pr(@T)’ 

A . 


a2 



Sodi urn 

0.0107(100) 

83.1 /AT^ 

Q.730 

1 .344 

0.975 

1 .012 

Potassium 

0.0085(100) 

73.5/AT 

0.511 

3.050 

0.975 

1 .011 

Bismuth 

0.0171(300) 


1 .74 

0.715 

1 .034 

1.20 

Lead 

0.0206(400) 


0.437 

2.20 

0.965 

1 .038 

Tin 

0.016(240) 

93/AT 

0.476 

2.32 

0.975 

1 .13 

Zinc 

0.0275(450) 

262/AT 

0.502 

2.41 

0.935 

1 .29 

Water 

13.6(0) 

1 59/AT 

0.255 

9.41 

1 .090 

2.008 


Notes: 1. Temperatures measured in °C. 


2 . Properties determined at Tm*for sodium and potassium and at 
25°C for the solid phase for other materials. 
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Table III 


EFFECT OF A ON CRITICAL MARANGONI NUMBERS 
FOR STEADY STATE PROCESS 

B=2.0, A=1.5, K=0.7 


A 

Merit 

ct (Merit) 

2 

142.6 

2.3 

20 

150.0 

2.4 

200 

150.6 

2.4 



TABLE IV 


Effect of A on Critical Marangoni Number 
for the Constant Rate Process 


A. Effect of Change of L/Cp^. 

Pr^O.Ol; K=0.7; Q=2.0; P=1.0; C=1.0; B=1 .0 


A 

R 

A^crit 

a(Afcrit) 

20 

14.0 

3.26x10^ 

2.1 

35 

8.2 

6.40x10^ 

2.1 

70 

4.2 

1 .61x10^ 

2.1 

105 

2.8 

2.76x10^ 

2.1 

175 

1 .7 

5.27x10^ 

2.2 


B. Effect of Change of • 

Pr=0.01; K=0.7; P=1 .0; C=1 .0; B=1 .0 


A 

_a_ 

R 

^crit 

^ C^crit) 

17.5 

0.5 

8.4 

6.22x10^ 

2.1 

35 

1 .0 

5.6 

1 .08x10^ 

2.1 

70 

2.0 

4.2 

1 .61x10^ 

2.1 

105 

3.0 

3.7 

1 .90x10^ 

2.1 

175 

5.0 

3.1 

2.18x10^ 

2.1 
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TABLE V 


Effect of Pr on Stability for R=1 .0 
Q=1.0; B=1.0; K=1 .0; P=1 .0; C=1.0 


Pr 

A 

^crit 

^crit 

a (A^crit) 

0.002 

1000 

4.96x10^ 

99.2 

2.2 

0.01 

200 

9.96x10^ 

99.6 

2.2 

0.1 

20 

1 .05x10^ 

101 .8 

2.2 

0.3 

6.0 

3.59x10^ 

107.7 

2.2 

0.6 

2.7 

2 

2.00x10 

120.0 

2.1 

0.9 

1 .5 

1 .45x10^ 

130.5 

1 .9 
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APPENDIX I 


Matching Conditions 

It is necessary to determine the conditions at the liquid-solid 
interface for arbitrary movement of that interface. Once the general 
relations have been derived the conditions will be linearized by 
assuming that the interface motion consists of a mean motion that 
is one-dimensional plus a small perturbation to the mean motion. 

The interface matching conditions v/ill be written for a fixed labora- 
tory frame of reference. 

General interface conditions, which are expressions of the con- 
servation of mass, momentum, and energy, are stated for a multi - 
component system by Williams (1965). These expressions can be special- 
ized for a control volume consisting of a thin slab that surrounds 
the interface. Assuming one-dimensional mean motion the slab control 
volume in a fixed coordinate system is shown in the figure below. 



Li qui d 


Sol i d 


S(t) - Interface 



Assume that there are no sources of mass at the surface, that the 

viscous stresses can be ignored, and that no chemical reactions occur. 

Write the velocities as follows: 

_» _fc. dS 

“s = dt = V 

, is the velocity of the control volume and k is a unit vector 

dt 

in the z direction. 

The continuity equation, for this control volume [Williams, p. 14], 


can be written 


Vdt^ ^ “ Ps dt ; 


Pa(v -w„) = p v 


The energy equation is 
.. 2 




where 


q = -k 


e = specific internal energy 
The momentum equation becomes 

PA^V''p) - Pd = ° 

Using this result and the continuity equation, the energy equation 
can be written as 

k; dT^ = Ps(®d-es> f 'Pd ^ '''p'”d^ 

d z j — dt 2 

dz 

But L=e„-e and the last term can be dropped because it is several 

36 S 

orders of magnitude smaller than the other terms. Thus, we get the 
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following equation that was used in the solution of the mean state for 
both problems. 


k dTs 
® ' 32 ' 



■3z“ 


PsLv 


P 


For the perturbation problems it was convenient to transform to a 
coordinate system fixed to the phase surface, S. If z denotes the 



This is eq. (16) that was used in the perturbation analysis. 
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Appendix II 

RESULTS FOR STEADY STATE PROBLEM 


A 

B 

Y 

A 

K 

M crit 

a(M crit) 

1000 

0.0 

0.2 

1 .5 

0.7 

79.603 

- 


0.5 

2.0 

0.2 



1 .5 



0.7 



145.7 


5.0 


150.2 


25.0 


150.6 


20.0 

0.0 

1 .2 

1 .5 

0.7 

79.0 

2.0 


2.0 


150.0 

wBsm 


4.0 


216.9 

2.5 


6.0 


282.5 

2.6 

2.0 

2.0 

1 .2 

1 .5 

0.7 

142.6 

2.3 

200.0 

2.0 

1 .2 

1 .5 

0.7 

150.6 

2.4 
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Appendix III 

RESULTS FOR THE CONSTANT RATE PROCESS 



R greater than 10 are in error because the linear equation was used for their determination 



- 2 - 


Pr 

A 

Q 

K 

P 

C 

B 

R 

Merit 

oc (Merit) 

0.002 

1000 

1 .0 

0.7 

1 .0 

1 .0 

1 .0 

1 .0 

4 . 96 x 10 ^ 

2.2 

0.01 

17.5 

0.50 

0.7 

0.97 

1 .0 

1 .0 

8.2 

6 . 60 x 10 ^ 

2.1 



0.50 


1.00 



8.4 

1 

6 . 28 x 10 ^ 

2.1 

0.01 

20.0 

2.0 

0.7 

1 .0 

1 .0 

1 .0 

( 

14.0 

3 . 29 x 1 0 ^ 

2.1 



3.0 





18.1 

2.44 

2.15 



4.0 





21 .9 

1 .96 

2.2 



5.0 





25.5 

1 .66 

2.2 

0.01 

35 

1 .0 

0.7 

1 .0 

1 .0 

1 .0 

5.6 

1 . 08 x 10 ^ 

2.1 



2.0 





8.2 

6 . 4 x 10 ^ 

2.1 

0.01 

70 

o 

CM 

0.7 

1 .0 

1 .0 

1 .0 

4.2 

1 . 61 x 10 ^ 

2.1 

0.01 

105 

2.0 

0.7 

1.0 

1 .0 

1 .0 

2.8 

2 . 76 x 10 ^ 

2.1 



3.0 





3.7 

1 .90 

2.1 

0.01 

175 

0.60 

0.7 

i,.o 

1 .0 

1 .0 

0.9 

1 . 10 x 10 ^ 

2.2 



0.65 





0.9 

1 .06 

2.2 



0.70 





1 .0 

1 .03 

2.2 



0.90 





1 .1 

9 . 05 x 10 ^ 

2.2 



1 .0 





1.1 

8.53 

2.2 



1 .1 





1 .2 

8.06 

2.2 



1 .2 





1 .2 

7.64 

2.2 



1 .5 





1 .4 

6.57 

2.2 1 



2.0 





1 .7 

5.27 

2.2 



2.5 





2.0 

4.36 

2.1 



3.0 





2.3 

3.69 

2.1 



3.5 




i 

2.5 

3.17 

2.1 1 



4.0 




1 

i 

2.8 

2.77 

2.1 ; 



4.5 




1 

1 

3,1 

2.45 

2.1 ; 



5.0 





3.4 

2.18 

2.1 j 

0.01 

175 

0.5 

0.7 

0.97 

1 .0 

1 .0 

0.8 

1 . 23 x 10 ^ 

i 

2.2 i 



1 .0 





1 .1 

8 . 89 x 10 ^ 

2.2 I 


i 

2.0 





1 .6 

5.51 

2.2 i 
1 


1 

1 

3.0 





2.2 

3.87 

2.1 ! 


1 

i 

! 

4.0 





2.7 

2.91 

2.1 i 


i 

1 

5.0 





3.3 

2.30 

2.1 j 
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-3 


Pr 

A 

Q 

K 

P 

C 

B 

R 

MO ri t 
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